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ABSTRACT 

We investigate the use of wavelet transforms in detecting and characterising non-Gaussian 
structure in maps of the cosmic microwave background (CMB). We apply the method to sim- 
ulated maps of the Kaiser-Stebbins effect due to cosmic strings onto which Gaussian signals 
of varying amplitudes are superposed. We find the method significantly outperforms standard 
techniques based on measuring the moments of the pixel temperature distribution. We also 
compare the results with those obtained using techniques based on Minkowski functionals, 
and we again find the wavelet method to be superior. In particular, using the wavelet tech- 
nique, we find that it is possible to detect non-Gaussianity even in the presence of a super- 
posed Gaussian signal with five times the rms amplitude of the original cosmic string map. 
We also find that the wavelet technique is useful in characterising the angular scales at which 
the non-Gaussian signal occurs. 

Key words: methods: data analysis - techniques: image processing - cosmic microwave 
background. 



1 INTRODUCTION 

The cosmic microwave background (CMB) is widely considered 
as one of the most important experimental tools in investigating 
the formation of structure in the Universe. In particular, by ob- 
serving the temperature fluctuations in the CMB, we hope to dis- 
tinguish between two competing theoretical paradigms. The first 
is the standard inflationary model, incorporating cold dark matter, 
for which the distribution of temperature fluctuations in the CMB 
should be Gaussian. The second class of theories invoke the for- 
mation of topological defects such as cosmic strings, monopoles 
or textures. In these theories the CMB temperature fluctuations are 
expected to be non-Gaussian, possessing steep gradients or well- 
defined 'hot spots' of emission (Bouchet, Bennett & Stebbins 1988; 
Turok 1996). Thus the detection (or otherwise) of a non-Gaussian 
signal in the CMB is an important means of discriminating between 
these two classes of theory. 

Even in defect theories, however, it is possible that the non- 
Gaussian CMB signal may be obscured by an additional Gaussian 
component of primordial CMB fluctuations, by foreground emis- 
sion from the Galaxy, or by instrumental noise. Moreover, CMB 
fluctuations arising from realistic defect cosmologies may well 
contain a Gaussian component that has a similar power spectrum 
to that of the non-Gaussian component (Magueijo & Lewin 1997), 
thereby making the task of detecting the non-Gaussian signal even 
more difficult. It is therefore important to investigate the regimes 
in which underlying non-Gaussianity can be detected when addi- 
tional Gaussian signals are present. The most common technique 
is to define a statistic with properties that are easily calculated for 
Gaussian processes and then apply this statistic to the (possibly) 



non-Gaussian CMB map. In this way it is often possible to deter- 
mine the probability that a given test map is drawn from an under- 
lying Gaussian ensemble. Several approaches to this problem have 
already been proposed. Perhaps the most straightforward of these is 
the measurement of the skewness and kurtosis of the distribution of 
pixel temperatures in the map (Scaramella & Vittorio 1991). More 
elaborate methods include the calculation of the 3-point correlation 
function of the temperature distribution (Kogut et al. 1996) and the 
investigation of the statistics of maxima and minima in the map 
(Bond & Efstathiou 1987; Vittorio & luskiewicz 1987). Methods 
based on the topological properties of the two-dimensional CMB 
temperature distribution have also been proposed (Coles 1988; Gott 
et al. 1990). The most recent of these techniques is based on the 
calculation of the three Minkowski functionals of the temperature 
distribution (see Section^)). 

In this paper, we investigate the use of a new method for 
detecting non-Gaussianity which is based on the idea of express- 
ing the CMB temperature distribution in terms of a set of two- 
dimensional wavelet basis functions (Fang & Pando 1996; Ferreira, 
Magueijo & Silk 1997). The wavelet transform provides a natu- 
ral decomposition of the image into structure on different scales 
and, by analysing the distribution of wavelet coefficients on a given 
scale, it is possible to avoid the usual restrictions imposed by the 
central limit theorem, so that non-Gaussian signal is more easily 
detected. We apply this technique to simulated maps of the CMB 
that contain both a non-Gaussian contribution, due to the Kaiser- 
Stebbins effect from cosmic strings (Kaiser & Stebbins 1984), and 
also a Gaussian component with an identical power spectrum to 
the strings map. By varying the relative amplitudes of the two con- 
tributions, we determine to what extent the non-Gaussian signal 
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can be obscured while still remaining detectable by the wavelet 
method. We compare our results with those obtained using either 
Minkowski functionals or straightforward calculation of the mo- 
ments of the pixel temperature distribution. 



2 THE DISCRETE WAVELET TRANSFORM 

The discrete wavelet transform (DWT) has been discussed exten- 
sively elsewhere (e.g. Daubechies 1992). In particular, a charac- 
teristically clear introduction is given by Press et al. (1994). We 
therefore give only a brief description of its properties. 



2.1 The one-dimensional DWT 

It is best initially to discuss the DWT in just one dimension. In 
order to develop an intuitive understanding, let us begin by consid- 
ering a periodic function f(x) of period L that we wish to expand 
in terms of a wavelet basis. The construction of the wavelet basis 
begins with the specification of the mother (or analysing) wavelet 
\|/(jc), together with another related function §(x) called the father 
wavelet (or scaling function). In order to construct a wavelet basis 
that is discrete, orthogonal and compactly-supported, the functions 
4»(jc) and y(x) must together obey several highly-restrictive math- 
ematical relations first derived by Daubechies (1992). We will not 
list these relations here, but merely note that the most straightfor- 
ward requirements are that 
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Mother and father functions with compact support are usually de- 
fined on the interval [0, 1], and the wavelet basis is then constructed 
from dilations and translations of §(x) and as 



h,i( x ) 



1/2 



ty(2 J x/L-I), 



1/2 



y(2 J x/L-l), 
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where j and I are integers. The index j labels the scale size of the 
wavelet, whereas / labels the position of the wavelet at this scale. 
The set {tyo.i,\fjj} with ^ j < °° and — «j < / < oo forms a com- 
plete, orthonormal basis in the space of functions of period L. The 
orthogonality of these functions is formally given by 

§jj(x)§j.l>(x)dx = 8«', 

(*)¥// = 5 //'' 8 H'> 

¥/,/(■*)¥// (*)d* = Sj/S//'- 

Using these orthogonality relationships, we may then write the 
function fix) as 

/(*)= E a 0,/<l>0./M+ E E bjjVjjix), (5) 
/=-» ;=0/=-oo 

where the wavelet coefficients a^ i and bj ; are given by 

ao,l = J f{x)$ 0t i(x)dx, (6) 



J_J{x)^jj{x)6x. 



(7) 



So far we have assumed our function f(x) is defined for all 
values of x. In many applications, however, the function is in fact 
digitised (sampled) at N = 2 J equally-spaced points Xj, In this case, 
by analogy with the Discrete Fourier Transform (DFT) of a sam- 
pled function, (pj) becomes 



y-i2'-i 

/(*;) =«o,o<t»o,o(^) + E E b j.Wj,i( x i. 

7=0 1=0 



(8) 



and the integrals for the wavelets coefficents in (g) & (|7|) are re- 
placed by the appropriate summations. If the mean of the function 
samples, /(x,), is zero (as will be case in this paper), then ao,o = 
and the function can be described entirely in terms of the wavelets 
V|/ ; - / . As the scale index j increases from to J — 1 , the wavelets rep- 
resent the structure of the function on increasingly smaller scales, 
with each scale a factor of 2 finer than the previous one. The index 
I (which runs from to 2 J — 1) denotes the position of the wavelet 
\\fjl within the jth scale level. 

If we collect the function samples /,- = /(jc;) in a column vec- 
tor f (whose length N must an integer power of 2), then (like the 
DFT) the DWT is a linear operation that transforms f into another 
vector f of the same length, which contains the wavelet coefficients 
of the (digitised) function. The action of the DWT can therefore be 
considered as the multiplication of the original vector by the N xN 
wavelet matrix W to give 



f = Wf. 



(9) 



Again like the DFT, the wavelet matrix W is orthogonal, so the 
inverse transformation can be performed straightforwardly using 
the transpose of W. Thus both the DFT and DWT can be considered 
as rotations from the original orthonormal basis vectors e, in signal 
space to some new orthonormal basis e; (i = l,...,N), with the 
transformed vector f containing the coefficients in this new basis. 

The original basis vectors e; have unity as the ith element and 
the remaining elements equal to zero, and hence correspond to the 
'pixels' in the original vector f . Therefore the original basis is the 
most localised basis possible in real space. For the DFT, the new 
basis vectors e,- are (digitised) complex exponentials and represent 
the opposite extreme, since they are completely non-local in real 
space but localised in frequency space. For the DWT, the new ba- 
sis vectors are the wavelets, which enjoy the characteristic property 
of being reasonably localised both in real space and in frequency 
space, thus occupying an intermediate position between the origi- 
nal 'delta-function' basis and the Fourier basis of complex expo- 
nentials. Indeed, it is the simultaneous localisation of wavelets in 
both spaces that makes the DWT such a useful tool for analysing 
data in wide range of applications. 

Since the wavelet basis e,- consists of (digitised) dilations and 
translations of the mother and father wavelets vy(x) and ty(x), we 
can obtain a different wavelet basis for each pair of these func- 
tions that obey the Daubechies conditions. Thus there exists an in- 
finite number of possible wavelet transforms, with different wavelet 
bases making different trade-offs between how compactly they are 
localised in either real space or frequency space. In this paper, we 
will be concerned only with wavelet bases for which the wavelets 
are continuous and have compact support in either space (it is im- 
possible for a function to have compact support in both spaces). 

The implementation of the DWT is based on a pyramidal al- 
gortihm (Daubechies 1992; Press et al. 1994), and its speed is com- 
parable to a Fast Fourier Transform. The calculation is arranged so 
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Figure 1. The domains of the matrix T containing the coefficients of the 
two-dimensional wavelet transform. The italic numbers show the value of 
k = 2 Ji + 2 ]2 in each domain; see text for details. 



that the transformed signal vector f = Wf , containing the wavelet 
coefficients, takes the form 

f = («o,o , ^o,o ,h,o,h,i, bi.o ,b2,i, ha , ^2,3 ,•■•)' • 

Thus the wavelet coefficients are placed in order of increasing scale 
index j and, within a given scale, the 2 J coefficients are ordered 
with increasing position index /. 



2.2 The two-dimensional DWT 

The extension of the DWT to two-dimensional signals (or images) 
is straightforward and is usually performed by taking simple tensor 
products of the one-dimensional wavelet basis to produce a two- 
dimensional wavelet basis given by 

<h,o;W2( x >>') = h,h( x )^o,h(y) 
ZjuWukfay) = VjiaM^oaCv) 
$oj 2 ;h,h(x,y) = <hA(jc)¥>2,feCy) 

If the two-dimensionalpixelised image has dimensions 2 J ' x 2 Jl 
then by analogy with (|8j) we have ^ j\ ^ J\ — 1 and ^ l\ ^ 
2 J ' — 1, and similarly for ji and 1%. 

In this paper, will we be concerned with two-dimensional pix- 
elised maps of temperature fluctuations in the CMB with dimen- 
sions 256 x 256 (i.e. J\ = J2 = 8). If we denote such an image by 
the matrix T, then the matrix of wavelet coefficients is given by 

f = WTW\ (10) 

where W is the wavelet matrix introduced in (^) and W' is its trans- 
pose. The structure of the matrix T containing the wavelet coef- 
ficients is shown in Fig. [j], where the pixel numbers are plotted 
on a logarithmic scale. We see that the matrix is partitioned into 
separate domains according to the scale indices j\ and j'2 in the 
horizontal and vertical directions respectively. By analogy with the 
one-dimensional case, as j\ increases the wavelets represent the 
horizontal structure in the image on increasingly smaller scales. 



Similarly, as ji increases the wavelets represent the increasingly 
fine scale vertical structure in the image. Thus domains that lie in 
the leading diagonal in Fig.|j| (i.e. with j\ = jj) contain coefficients 
of two-dimensional wavelets that represent the image at the same 
scale in the horizontal and vertical directions, whereas domains 
with j\ ^ J2 contain coefficients of two-dimensional wavelets de- 
scribing the image on different scales in the two directions. 

Finally, we note that for the two-dimensional DFT it is com- 
mon to denote the distance from the origin in Fourier space by 
k, which serves as a measure of inverse scale-length. In a two- 
dimensional DWT, it is also useful to define a similar quantity, 
but in this case k is an integer variable given by k = 2 JI + 2 n 
which takes certain prescribed values between k m [ a = 2 and k m!a = 
2 J, ~ l +2 Jl ~ l . In our application J\ = J2 = 8 and so k max = 256. 
By definition, the value of k is constant within each of the domains 
shown in Fig|l| The values of k in each domain are shown as the 
italic numbers in Fig. |l| We note that k takes the same value for do- 
mains lying symmetrically on either side of the leading diagonal. 



3 DETECTING NON-GAUSSIANITY WITH THE DWT 

As mentioned above, let us denote our two-dimensional pixelised 
test map of temperature fluctuations in the CMB by the matrix T. If 
the test map is not drawn from a Gaussian ensemble, then we might 
hope to measure this non-Gaussianity directly from the histogram 
of the pixel temperatures, i.e. the histogram of the elements of the 
matrix T. Indeed, Scaramella & Vittorio (1991) use the skewness 
and kurtosis of the pixel temperature distribution as statistics for 
detecting non-Gaussianity in CMB images. 

As noted by Fang & Pando (1996), however, methods based 
on statistics of the pixel temperature distribution can be unreliable 
in detecting non-Gaussianity, as a result of the central limit theo- 
rem. For example, if an image consists of numerous non-Gaussian 
processes on different scale lengths, then from the central limit 
theorem we would expect the probability distribution of the pixel 
temperatures to tend to a Gaussian as the number of processes in- 
creases. 

In this paper, therefore, rather than analysing the distribution 
of the elements of T, we concentrate instead on the statistics of 
the wavelet coefficients, contained in the matrix T. The wavelet 
transform provides a natural decomposition of the image into struc- 
ture on different scales, and within each scale the wavelet basis 
functions are well-localised. Thus, by analysing the distribution of 
wavelet coefficients on a given scale, it is possible to avoid the re- 
striction of the central limit theorem and any non-Gaussian signal 
on this scale is far more pronounced. 



3.1 Statistics of wavelet coefficients 

Referring to Fig. [lj it is useful to consider together all the wavelet 
coefficients that lie in domains sharing the same value of k. It is 
straightforward to show (and is intuitively reasonable) that if the 
input image T is (cyclically) translated or rotated, then wavelet co- 
efficients sharing the same lvalue are merely rearranged within 
the relevant domains. Thus any statistic based on the distribution 
of all wavelet coefficients sharing the same fc-value is invariant to 
translations and rotations of the input image. 

In order to describe the statisics of the wavelet coefficients 
of a given test image, we could (for each value of k separately) 
calculate estimators of the moments v r (r = 1,2,3,.. .), or central 
moments fi r , of the distributions of its wavelet coefficients (see the 
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Appendix). In this way, we obtain moment spectra v r (k) (or central 
moment spectra fi r (k)) for r = 1,2,3,... which describe the (scale 
dependent) statistics of the wavelet coefficients. As discussed by 
Ferreira et al. (1997), however, a better approach is to describe the 
distribution of wavelet coefficients at each A>value in terms of its 
cumulants. Using unbiassed estimators of the cumulants based on 
fc-statistics (see the Appendix), we thus calculate cumulant spectra 
k r (k) (r = 1,2,3, ...) of the wavelet coefficients of the test image. 

For a given test image T, the cumulant spectra of the corre- 
sponding wavelet coefficients have some useful properties. Firstly, 
let us suppose that T is drawn from a Gaussian ensemble. Since the 
wavelet transform can be considered merely an orthogonal rotation 
of the basis vectors in image space, it is straightforward to show 
that the wavelet coefficients at each value of k are also drawn from 
a Gaussian distribution. Thus, for Gaussian test images, the expec- 
tation value of the estimated cumulant spectra k r (k) of the wavelet 
coefficients is zero for r > 2 (see the Appendix). Secondly, let us 
consider a test image T = Tj + T2 that is the sum of two indepen- 
dent processes T; and T; . Since the DWT is a linear operation, we 
thus have T = T~i + T2. Using the additive property of cumulants 
discussed in the Appendix, we therefore find that the rth cumulant 
spectrum k r (k) of the wavelet coefficients of T is simply the sum 
of the corresponding cumulant spectra of the wavelet coefficients 
of T~i and T2 taken separately. Of course, this argument can be 
extended to an arbitrary number of independent processes. 

The cumulant spectra k,-(k) (r = 1,2,3,...) form our basic 
test for non-Gaussianity. The procedure is as follows. We begin 
with some (possibly non-Gaussian) test image T and perform a 
wavelet transform to obtain the corresponding matrix of wavelet 
coefficients T given by (|Io|). We then calculate the cumulant spec- 
tra K r (k) of these wavelet coefficients. Returning to the original 
test image T, we then create a large number (5000 say) of equiva- 
lent Gaussian realisations (EGR). Each EGR is created by Fourier 
transforming the original test image, randomising the phases of the 
complex Fourier coefficients, and inverse Fourier transforming the 
result. In the randomisation step, the phase of each Fourier coeffi- 
cient is drawn at random from a uniform distribution in the range 
[0, 2tz] . Thus each EGR has exactly the same power spectrum as the 
original test image, but the pixel temperature distribution is drawn 
from a Gaussian ensemble. Each EGR is then analysed in exactly 
the same way as the original test image in order to obtain the cumu- 
lant spectra of its wavelet coefficients. Thus, after analysing all the 
EGRs, we obtain (for each value of k) approximate probability dis- 
tributions of the cumulant estimators k r (k) for a Gaussian process 
with the same power spectrum as the original test map. By com- 
paring these probability distributions with the cumuluant spectra 
of the original image, we thus obtain an estimate of the probabil- 
ity that the original image was drawn from a Gaussian ensemble. 
For 256 x 256 images considered here, the entire analysis requires 
about 20 minutes CPU time on a Sparc Ultra workstation. 

We note that, in two-dimensional applications it is usual to 
consider only the domains in Fig. [I] with 71,72 ^ !• The main 
reason is that the domains with j\ = or 72 = contain a mix- 
ture of coefficients of the three different two-dimensional wavelet 
bases l^o;/, {x,y), £7, ,0;h ,/ 2 (x,y) and ^ j 2; /, ,l 2 (x,y). Domains with 
./li 72 ^ 1. however, contain only the coefficients fcj]./ 2 ;/,,/ 2 of the 
basis x Vj l ,j 2 ;h,h( x 'y) anc ^' as a resu lt, statistics based on the distri- 
bution of these wavelet coefficients are generally more informative. 
Therefore, for the remainder of this paper, we will restrict our at- 
tention to domains with 71 ,72 ^ 1- 



3.2 Choosing the wavelet basis 

Our only remaining task is to choose the wavelet basis in which we 
wish to expand our CMB map. As mentioned above, in this paper 
we restrict ourselves to discrete, orthogonal, compactly-supported 
wavelet bases. Nevertheless, there are still an infinite number of 
such bases and so we will concentrate on only a small subset of 
commonly-used wavelet bases. Our wavelet 'library' consists of 
the following one-dimensional wavelet bases: Haar; Daubechies 
4,6,12,20; Coiflet 2,3; Symmlet 6,8. Several of these wavelets are 
plotted and discussed by Graps (1995). We note that this library 
is in no way intended to be an exhaustive set of useful wavelets, 
but represents a fairly typical range of bases which could be em- 
ployed. As mentioned above, in each case, the corresponding two- 
dimensional wavelet basis is formed by taking the tensor product 
of the one-dimensional basis. 

Since each wavelet basis provides a different decompostion of 
the CMB test image, the resulting cumulant spectra will be differ- 
ent in each case. Indeed, depending on the structure present in the 
test image, some wavelet bases will be more powerful at detecting 
non-Gaussianity than others. Since the analysis discussed above re- 
quires the wavelet transform of a large number of EGR to be calcu- 
lated, it would clearly be very time-consuming to perform this for 
every available wavelet basis. Therefore, for a given test image, we 
wish to find a quick method by which we can determine the wavelet 
basis that is most sensitive to any non-Gaussian structure present. 

We begin by calculating the wavelet coefficients by,.;,;/,,/, of 
the single test map in each of the available bases, which requires 
only a couple of seconds of CPU time on a Sparc Ultra workstation. 
For each basis, we then define the 'normalised' wavelet coefficients 
at each value of k by 



where, in the denominator, the sum on the scale indices j\ and 72 
is for all values satisfying 2" + 2 J1 =k and the sum on position in- 
dices l\ and I2 corresponds to all wavelets coefficients lying within 
the relevant scale domains. Thus, for each value of k, the coeffi- 
cients Pj l j 2 ;hh nave values between zero and unity and the sum 
of all the coefficients is equal to unity. We then calculate (a multi- 
ple of) the entropy of the normalised wavelet coefficients for each 
value of k, 

5= ~In~/V ^ PjuKh,h^Pjuh;h,h 
h,hhh 

where again the sum extends over all wavelet coefficients sharing 
the fc-value under consideration and N is the total number of these 
coefficients. We also note that lim^o plnp = 0. The value of the 
sum lies between and — In N, and so we have normalised 5 so that 
its value always lies between zero and unity. The entropy 5 takes 
the value zero if one of the coefficients Pj l .j 2 -j 2 ,i 2 equals unity and 
the rest are zero. At the other extreme, S equals unity when all the 
coefficients are equal. 

Now, if the original CMB map is Gaussian, then, as mentioned 
above, the distribution of wavelet coefficients at each k- value is also 
Gaussian. Furthermore, for a Gaussian CMB map with structure on 
a wide range of scales, we would expect that, at each scale k, all 
the wavelet basis functions would be required, with roughly equal 
amplitudes, in order to represent the image on that scale. Hence 
the entropy of the normalised wavelet coefficients should be close 
to its maximum value of unity. For a map with pronounced non- 
Gaussian features, however, this is not necessarily the case. If there 
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Figure 2. (a) A realisation of temperature fluctuations in the CMB due to 
the Kaiser-Stebbins effect from cosmic strings. The map is pixelised on a 
256 x 256 grid with a pixel size of 1.5 arcmin. For convenience the map 
has been normalised so that its variance equals unity, (b) A histogram of the 
pixel temperature distribution. 

exists a wavelet basis that resembles the non-Gaussian features in 
the map on some scale, then we might expect the original image to 
be well represented by fewer basis functions. In this case, many of 
the wavelet coefficients may be close to zero, with just a few large 
coefficients corresponding to the basis functions that best resemble 
the main features in the image. We would therefore expect a lower 
value for the entropy of the normalised wavelet coefficients. Thus, 
for each wavelet basis in our library, we calculate the entropy of 
the wavelet coefficients of the test image at each value of k. The 
optimal basis is then taken as that having the lowest entropy for any 
fc-value. For each of the test maps discussed in the next Section, the 
Coiflet 2 wavelet basis was found to be optimal. In particular, the 
entropy of the wavelet coefficients was especially low for k = 66, 
which, as we shall see below, is the angular scale on which a non- 
Gaussian signal is detected with the highest significance for each 
test map. 



4 APPLICATION TO SIMULATED CMB MAPS 

In this Section, we apply the wavelet transform technique described 
above to the detection of non-Gaussianity in simulated CMB maps. 



Figure 3. (a) A realisation of a Gaussian image with a Dower spectrum 
identical to that of the cosmic strings map shown in Fig. EJa); see text for 
details, (b) A histogram of the pixel temperature distribution. 



Our basic non-Gaussian test map is shown in Fig. |2|(a) and is a 
realisation of CMB anisotropies due to the Kaiser-Stebbins effect 
from cosmic strings (Bouchet, Bennett & Stebbins 1988). The map 
is pixelised on a 256 x 256 grid with a pixel size of 1.5 arcmin. 
Thus the size of the simulated field is 6.4 x 6.4 deg 2 . The map has 
zero mean and has, for convenience, been normalised so that its 
variance equals unity. 

We see from Fig. ^(a) that the strings map is extremely non- 
Gaussian, possessing sharp temperature gradients and localised 
hot-spots. From the scale-bar, we also note that the maximum and 
minimum temperatures in the map are respectively 6.1 and —5.3 
times the rms value, again indicative of a highly non-Gaussian pro- 
cess. This non-Gaussianity is also apparent from the histogram of 
pixel temperatures shown in Fig. ^(b). 

In order to assess to what extent the non-Gaussian signal in 
the strings map can be obscured, but still remain detectable, we will 
add to it some multiple of the Gaussian test map (GTM) shown in 
Fig. ^(a). The GTM shown is an equivalent Gaussian realisation of 
the cosmic strings map shown in Fig. ^(a) and thus has an identi- 
cal power spectrum, but the pixel temperatures are drawn from a 
Gaussian distribution. Clearly, by construction, the GTM also has 
unit variance. 

We could, of course, use a GTM that is not an equivalent Gaus- 
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sian realisation of the strings map. For example, we could use a 
realisation of CMB fluctuations predicted from some inflationary 
model. Such a map would, however, in general possess a different 
power spectrum to the strings map and differences between the two 
could be discerned by measuring these different power spectra. By 
using a GTM that is an EGR of the strings map, we are consider- 
ing the case of optimal confusion between the non-Gaussian and 
Gaussian signals and this should provide the most stringent test 
of the methods for detecting non-Gaussianity. Moreover, as men- 
tioned in the Introduction, CMB fluctuations from realistic cosmic 
string cosmologies may well contain a Gaussian component that 
has a similar power spectrum to that of the non-Gaussian compo- 
nent. (Magueijo & Lewin 1997). 

The differences between the strings map and the GTM are 
clear to the eye. The GTM possesses no sharp features or extreme 
maxima or minima. Furthermore, from the scale-bar, we see that 
the maximum and minimum of the pixel temperature distribution 
lie at 3.6 and -3.8 times the rms level, as one might expect for a 
Gaussian image with this number of pixels. The histogram of the 
pixel temperatures in the GTM is shown in Fig. ^(b) and is clearly 
closer to a Gaussian in shape than the corresponding histogram for 
the strings map shown in Fig. [^b). 

Our test maps thus consist of the sum of the above strings 
map and GTM in various proportions (a: b). In order to make the 
simulations more realistic, before being analysed the test maps are 
first convolved with a 5-arcmin Gaussian beam and Gaussian pixel 
noise is added with an rms level equal to one-tenth that of the con- 
volved map. This degree of smoothing and noise level is typical 
of what might be achieved for CMB observations using the ESA 
Planck Surveyor satellite (Bersanelli et al. 1996), after foreground 
contaminants have been removed using the maximum-entropy sep- 
aration algorithm discussed by Hobson et al. (1998). The process 
of convolution smoothes out many of the sharp features visible in 
the underlying cosmic strings map in Fig. |^(a), and the addition of 
pixel noise obscures this structure still further. In each case, the re- 
sulting convolved, noisy test map is rescaled to have unit variance 
before being analysed for any non-Gaussian signal. 

As mentioned above, for a wide range of non-Gaussian-to- 
Gaussian proportions (a:b), the Coiflet 2 wavelet basis was found 
to be optimal for detecting the non-Gaussian signal due to the 
strings map in Fig. ^](a). Therefore, for the remainder of this pa- 
per, all wavelet decompositions will be in this basis. Furthermore, it 
was found, for our simulations, that the fourth cumulant spectrum, 
K4(k), was the most sensitive to the presence of any non-Gaussian 
signal, and so we will plot only this spectrum for each test map. 



4.1 The Gaussian test map 

We begin by assuming the CMB fluctuations to consist only of 
the Gaussian map shown in Fig. ^(a). Thus, formally, our initial 
test map consists of the sum of the cosmic strings map and the 
GTM in the proportions (0:1) . This is then convolved with a 
5-arcmin beam and pixel noise is added as discussed above. Ap- 
plying the wavelet non-Gaussianity test to this image provides a 
useful check that a detection of non-Gaussianity is not obtained 
for a purely Gaussian signal. The resulting fci(fc) cumulant spec- 
trum of the wavelet coefficients of the test map is plotted as the 
solid squares in Fig. Q As mentioned in Section p\l[ the same non- 
Gaussianity test is then also applied to 5000 equivalent Gaussian 
realisations (EGR) of the test map. In Fig. ^, the mean value of 
K4(fc) obtained from the 5000 EGR is shown by the open circles, 
which all lie close to zero, as expected for a Gaussian process. At 



Figure 4. The K4(k) spectrum (solid squares) of the wavelet coefficients 
of the test map consisting of the sum of the cosmic strings map shown in 
Fig. ^(a) and the Gaussian test map shown in Fig. ^(a) in the proportions 
(0 : 1), after convolution with a 5-arcmin Gaussian beam and the addition of 
Gaussian pixel noise with rms equal to one-tenth that of the convolved map. 
The open circles show the mean ict(ft) spectrum obtained from 5000 equiv- 
alent Gaussian realisations of the test map (see text for details) and the error 
bars show the 68, 95 and 99 per cent limits of the resulting K4 distribution at 
each value of k. For presentational purposes, the spectrum been normalised 
so that, at each value of k, the variance of the K4 distribution obtained from 
the 5000 EGR is equal to unity. The point at k = 1 is not calculated from the 
wavelet coefficients, but instead shows the values of K4 obtained directly 
from the pixel temperature distributions in the test map and the 5000 EGR. 



each value of k, the error bars show the 68, 95 and 99 per cent limits 
of the K4 distribution obtained from the 5000 EGR. For presenta- 
tional purposes, the cumulant spectrum has in fact been normalised 
so that, at each value of k, the variance of the K4 distribution ob- 
tained from the 5000 EGR is equal to unity. Thus the significance 
of any detection of non-Gaussianity can be read off directly from 
the scale on the vertical axis of the plot. As might be expected in 
this case, however, the fci(fc) spectrum of the test map lies well 
within the 2-sigma limits of the Gaussian ensemble at all values of 
k. ____ 

since we ignore domains of the 



As discussed in Section 3.1, 
wavelet transform with j\ or 72 equal to zero, we see from Fig. |l| 
that the iCt(fc) spectrum of the wavelet coefficients should start at 
k = 4. The point plotted at k = 1 in Fig. ^ is not in fact calcu- 
lated from the wavelet coefficients, but instead shows the values of 
£4 obtained directly from the pixel temperature distributions in the 
test map and the 5000 EGR. This point therefore represents a non- 
Gaussianity test similar to that proposed by Scaramella & Vittorio 
(1991), in which the skewness and kurtosis of the pixel tempera- 
ture distributions were calculated. The point is included to provide 
an immediate comparison between methods based on the statistics 
of the wavelet coefficients and those based directly on the statistics 
of the pixel temperature distribution. 



4.2 The cosmic strings map 

We may repeat the above analysis for the opposite extreme in which 
the CMB fluctuations are assumed to consist only of the cosmic 
strings contribution (so that the non-Gaussian:Gaussian propor- 
tions are formally (1 : 0)). As before, the map is then convolved 
with a 5-arcmin beam and pixel noise is added. The resulting fct(fc) 
spectrum is shown in Fig. pL We see from the figure that unambigu- 
ous detections of non-Gaussianity are obtained at numerous values 
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Figure 5. As in Fig. ^, but for the test map with non-Gaussian: Gaussian 
proportions (1 : 0). 



of k, the most spectacular being a 910-sigma detection at k = 66. 
As discussed in Section 3.2, the entropy of the wavelet coefficients 
for the Coiflet 2 basis was lowest in the domains with k — 66, and 
so we might have expected these domains to provide the most sig- 
nificant detection of non-Gaussianity. 

Indeed, since the detection of non-Gaussianity is so large, it is 
not possible to discern the error bars denoting the confidence limits 
of the K4 distributions for the 5000 EGR (although they are plotted). 
We merely note that for the k — 1 point (calculated directly from the 
pixel temperatures), the significance of the non-Gaussian detection 
is only 1.1-sigma. Thus, in this case we find that, by analysing the 
wavelet coefficients of the image in separate /^-domains, we obtain 
a vastly improved detection of non-Gaussianity as compared with 
analysing the pixel temperature distribution directly. 



4.3 Mixed cosmic strings/Gaussian maps 

Now we have analysed separately the cosmic strings map in 
Fig. ^a) and the GTM in Fig. ^|(a) (after convolution and the ad- 
dition of pixel noise), we may repeat the analysis for test maps 
consisting of the sum of the cosmic strings map and GTM in dif- 
ferent proportions (1 : b). It is clear that when b 1 the resulting 
K4(k) spectrum of the wavelet coefficients will tend to that shown 
in Fig. W. Similarly, if b <C 1 the K^{k) spectrum will lie close to 
that in Fig. ^|. In this Section, our aim is to find the value of b at 
which the detection of non-Gaussianity becomes marginal. We re- 
iterate that, in the following analysis, each test map is convolved 
and has noise added (as discussed above) before applying the non- 
Gaussianity test. 

In Fig. ^ we plot the i<4(£) spectra for test maps with b — 
1,2,3,4,5,6. We see from the figure that significant non-Gaussian 
detections are obtained for b ^ 4. We note that, in each case, the 
largest detection of non-Gaussianity occurs at k = 66, as expected 
from the entropy analysis discussed in Section 3.2 The significance 
of the largest detections for b = 1,2,3,4 is 495-, 88-, 24- and 9- 
sigma respectively. For b = 5, the detection of non-Gaussianity be- 
comes marginal. The largest signal again occurs at k = 66, but the 
significance is only 4-sigma. For b = 6, no detection is possible 
and, as expected, the lct(fc) spectrum closely resembles that shown 
in Fig. ^ for the case when no non-Gaussian signal is present. 

For all values of b, however, we note that for the k = 1 point, 
which is calculated directly from the pixel temperatures, the value 



of K4 for the test map lies within the 95 percent confidence limits 
of the of K4 distribution obtained from the 5000 EGR. Thus, we see 
the statistic based on the pixel temperature distribution is far less 
sensitive to the presence of a non-Gaussian signal than those based 
on the wavelet coefficients of the image. 

The performance of the wavelet non-Gaussianity test is en- 
couraging. In particular, we must remember that the parameter b 
represents the ratio of the rms values of the GTM and cosmic 
strings maps that constitute the test image. The ratio of the power 
in the Gaussian and non-Gaussian contributions is given by b 2 and 
so a marginal detection is still possible when the power due to the 
GTM is 25 times that due to the cosmic strings map. We note that, 
since the GTM possess an identical power spectrum to the cos- 
mic strings map, the ratio of the power in the Gaussian and non- 
Gaussian contributions is equal to b 2 on all angular scales. Thus, 
the detection of non-Gaussianity is not the result of the presence of 
excess power in the cosmic strings map over some range of scales. 

For illustration purposes, in Fig. |v|(a) we plot the b = 5 test 
map, after convolution and the addition of pixel noise. Although 
the wavelet non-Gaussianity test provided a marginal detection, as 
shown in Fig. [](e), we see that, at least by eye, the non-Gaussian 
signal from the cosmic strings map in Fig. Pla) is completely ob- 
scured by the signal due to the GTM in FigTpla). Moreover, con- 
volution with a 5-arcmin Gaussian beam and the addition of pixel 
noise have further diluted the non-Gaussian signal. The histogram 
is pixel temperatures is shown in Fig. f^(b) and is clearly close to 
Gaussian. 



4.4 The cosmic strings map with enhanced pixel noise 

As a final example, we consider again the case in which the CMB 
signal consists only of the cos mic strings map, convolved to 5- 
arcmin resolution (as in Section 4.2 ), but in which the level of in- 



strumental pixel noise is increased such that its rms value is equal 
to that of the convolved CMB map. This example is included for 
two reasons. Firstly, it has traditionally been pixel noise that pre- 
sented the greatest obstacle to the detection of non-Gaussianity in 
the CMB, and this is still the case for the 4-year COBE data. Al- 
though the observational strategy for our simulated data is very dif- 
ferent from that used for the COBE observations, the information 
content of the two data-sets is comparable. The 6.4 x 6.4 deg 2 map 
used here contains approximately 6000 5-arcmin beams, which is 
similar to the number of pixels in standard resolution COBE ob- 
servations, although this somewhat larger than the number of in- 
dependent COBE beams. The second reason for including this test 
image is to investigate the abilities of the wavelet technique when 
the main contaminant does not possess a power spectrum that is 
identical in shape to that of the non-Gaussian signal, but instead 
has a flat white-noise power spectrum. 

The resulting i<4(fc) spectrum for this test map is shown in 
Fig. |^. We see from the figure that the most significant detection of 
non-Gaussianity is 42-sigma at k = 34, accompanied by a 24-sigma 
detection at k = 66 and several 5-10-sigma detections at other k- 
values. It is interesting to note that the most significant detection 
is much smaller than the 495-sigma detection obtained in Fig. ^(a) 
and occurs at a lower value of k. For Fig. ^J(a), the main contam- 
inant also had an rms equal to that of the cosmic strings map, but 
possessed an identical power spectrum. It is clear therefore that the 
wavelet technique is less effective at detecting the non-Gaussian 
signal due to the cosmic strings map when the main contaminant 
has a flat white-noise power spectrum. This is in fact to be ex- 
pected, since the non-Gaussian signal in the strings map comes 
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Figure 6. As in Fig. ^[ but for test maps with non-Gaussian: Gaussian proportions (1 : b). 



predominantly from sharp temperature gradients, which occur on 
small scales. As compared to the cosmic strings map, a pixel noise 
contaminant with the same rms value has relatively more power 
on small scales, and is therefore more effective at obscuring the 
non-Gaussian signal. By comparing Figs ^(a) and ^ more closely, 
this behaviour becomes more apparent. We see that in both cases a 
~ 40-sigma detection is obtained at k = 34, but the more significant 
detections in Fig. ^a), which occur at higher ^-values, have been 
much reduced in Fig. ^. 



5 MINKOWSKI FUNCTIONALS ANALYSIS 

In the previous Section, we found that the wavelet non-Gaussianity 
test was far more sensitive to the presence of a non-Gaussian signal 
than statistics based simply on the pixel temperature distribution of 
the CMB map. In order to compare the wavelet approach with other 
methods for detecting non-Gaussianity, in this Section we reanalyse 
the above test maps using an approach based on the topology of the 
CMB fluctuations. 

The approach we adopt here requires the calculation of three 
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Figure 7. (a) The test map containing.^ the cosmic strings map shown in 
Fig. |^(a) and the EGR shown on Fig. g(a) in the proportions (1 : 5). The 
map has been convolved with a 5-arcmin Gaussian beam and Gaussian pixel 
noise has been added with an rms equal to one-tenth that of the convolved 
CMB map. For convenience the map has been normalised to that its variance 
equals unity, (b) A histogram of the pixel temperature distribution. 



fll 



Figure 8. As in Fig. ^, but for the test map with non-Gaussian: Gaussian 
proportions (1 : 0) to which Gaussian pixel noise has been added with an 
rms value equal to that of the convolved CMB map. 



Minkowski functionals for each test map. Minkowski function- 
als have been applied to CMB data in several forms previously. 
Schmalzing & Gorski (1998) derive the Minkowski functionals 
for a spherical sky and apply their results to the COBE four year 
data. Winitzki & Kosowsky (1998) give a detailed derivation of the 
expected form of the Minkowski functionals for a flat Euclidean 
space, taking into account boundary conditions, pixel noise and 
pixel size and shape. The reader is refered there for a more the- 
oretical overview. 

We opt here for an approach that is as simple as possible and 
use the three two-dimensional Minkowski functionals on square 
pixels in flat space. This approach is sufficient for the compari- 
son desired in this paper. The three Minkowski functionals are the 
surface area, perimeter and Euler characteristic of an excursion re- 
gion. The excursion region R(T) is taken as the region of the map 
above a certain threshold temperature T, which is varied between 
the minimum and maximum values in the test map. The Minkowski 
functionals are therefore functions of the threshold temperature T. 
In two dimensions, the three functionals completely describe the 
statistical properties of a distribution. In terms of integral geometry 
they are given by 



dS, 



R(T) 



dl. 



M (T) = / 

■I R 

Mi(T) = [ 

J3R(T) 

M 2 (T) = / -. 

where dR(T) is the boundary of the excursion region R(T) at the 
threshold temperature T . The differentials dS and dl denote respec- 
tively the elements of area and of length along the boundary, and r 
is the radius of curvature of the boundary. As an example, a circular 
disc of radius a would have Mq = Tta 2 , Mi = 2na and M 2 = 2ti. 

The surface area and perimeter of the excursion region in a 
pixelised map are trivial to calculate but the Euler characteristic re- 
quires a little explanation. It should be noted that, as defined here, 
the third Minkowski functional is equal to 2nG where G is the 
genus of the excursion region. The genus is defined as the num- 
ber of isolated holes (regions below the threshold) minus the num- 
ber of islands (regions above the threshold) within the map. Using 
the approach advocated by Melott et al. (1989), based upon the 
three dimensional genus calculation of Hamilton, Gott & Wein- 
berg (1986), we calculate the genus of the two-dimensional map 
by looking at the angle deficits of the vertices. Using this method it 
is also straightfoward to calculate the three Minkowski functionals 
simultaneously. 

This approach to genus has been implemented previously and 
applied to the COBE data by Smoot et al. (1994), Kogut et al. 
(1996) and Colley, Gott & Park (1996). Gott et al. (1990) also 
looked at the use of the perimeter of the excursion regions but did 
not apply this analysis to any data. Theoretical predictions from de- 
fect theories have also been tested using the genus algorithm (see 
e.g. Avelino et al. (1998)) but the full Minkowski functional set 
is yet to be applied to a full CMB data set to constrain any cos- 
mological parameters or defect theories. This has been carried out 
for various Large Scale Structure data sets (see e.g. Kerscher et al. 
1997) with promising results. 
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Figure 9. Minkowski functionals analysis of the test map with non-Gaussian: Gaussian proportions (1,0). Top panels: the three Minkowski functionals of the 
test map (dashed lines) calculated at 50 separate threshold temperatures between the minimum and maximum of the test map pixel temperature distribution. 
The solid lines show the mean Minkowski functionals calculated from 5000 EGR and the error bars denote the one-sigma limits of the probability distributions 
of the functionals at each temperature step. Bottom panels: the difference between the Minkowski functionals of the test map and mean Minkowski functionals 
of the 5000 EGR. At each temperature step, the plots have been normalised so that the one-sigma points of the probability distributions obtained from the 5000 
EGR are equal to unity. 



5.1 Application to simulated CMB maps 

The three Minkowski functionals were calculated for the simulated 
CMB maps analysed in the previous Section. As discussed above, 
these test maps consist of the strings maps in Fig. ^(a) and the EGR 
in Fig ^(a) in proportions (1 : b), where b can vary. Before being 
analysed, the test maps are first convolved with a 5-arcmin Gaus- 
sian beam and Gaussian pixel noise is added with an rms equal to 
one-tenth the rms of the convolved map. In each case the tempera- 
ture range of the test map is divided into 50 steps and, at each tem- 
perature step, the three Minkowski functionals of the correspond- 
ing excursion region R(T) are calculated. Following the procedure 
adopted in the wavelets non-Gaussianity test, in each case the three 
Minkowski functional are also calculated in an analogous manner 
for 5000 EGR so that approximate probability distributions are ob- 
tained for the values of the functionals at each temperature step. 

The results for the cases b = 0, 1 and 2 are shown in Figures 
|ll| In the top three panels of each figure, we plot the Minkowski 
functionals of the test map (dashed lines) and the mean Minkowski 
functional for the 5000 EGR (solid lines). The error bars in each 
panel correspond to the one-sigma limits of the probability distri- 
butions of the values of the functionals at each temperature step, as 
obtained from the 5000 EGR. In order to estimate the significance 
of any non-Gaussian detection, in the bottom three panels of each 
figure we plot the difference between the Minkowski functionals of 
the test map and mean Minkowski functionals of the 5000 EGR. 
At each temperature step, the plots have been normalised so that 



the one-sigma points of the probability distributions obtained from 
the 5000 EGR are equal to unity. Thus, the significance of any non- 
Gaussian detection may be read-off directly from the vertical axis 
of each plot. 

In the b = case (for which the test image is simply the con- 
volved strings map with noise added), we see from the bottom three 
panels in Fig. [] that significant detections of non-Gaussianity are 
obtained for all three Minkowski functionals. For low and high val- 
ues of T, the corresponding excursions region of the test map are 
significantly different from those in a typical EGR and all three 
functionals detect some non-Gaussianity. Indeed the Mo, M\ and 
Mi functionals have detections of 14, 20 and 24-sigma respectively 
at high positive threshold temperatures. For the M\ functional, we 
also see that, when the threshold temperature T is close to zero 
(the mean of the maps), the excursion region R(T) of the test map 
has a significantly smaller perimeter than for a typical EGR and a 
6-sigma detection of non-Gaussianity is obtained. 

This behaviour is easily understood by examining the structure 
present in the cosmic strings map shown in Fig. ^. The strings map 
possesses a great deal of structure for threshold temperatures near 
the minimum and maximum of the pixel temperature distribution. 
Even after convolution with the 5-arcmin beam, there still exist sev- 
eral sharply defined hot-spots and cold-spots that lie outside the 
typical 3-sigma range of the corresponding EGRs, and hence lead 
to large detections of non-Gaussianity. As noted earlier, the range of 
values in the (unconvolved) strings map is 6.1 to -5.3, and so at its 
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Figure 12. As in Fig. Q but for the test map with non-Gaussian: Gaussian proportions ( 1 : 0) to which Gaussian pixel noise has been added with an rms value 
equal to that of the convolved CMB map. 



extremes the temperature distribution is slightly skewed to positive 
values, which leads to more pronounced detections at high tem- 
peratures. Thus, as might be expected, the Minkowski functional 
analysis is powerful at detecting non-Gaussianity due to points in 
the strings map temperature distribution that lie outside the corre- 
sponding distribution of a typical EGR. For threshold temperatures 
near zero, however, the strings map is quite structureless, possess- 
ing large plane areas at constant temperature. Thus, the perimeter of 
a correpsonding excursion region will be quite small in comparison 
to an EGR and leads to the 6-sigma detection in the Mi functional. 

From Figs [l(] & [llj however, we see that the results for 
b = 1 and b = 2 are somewhat disappointing. None of the 
Minkowski functionals have produced a significant detection of 
non-Gaussianity. Very few points lie outside the one-sigma limits 
obtained from the 5000 EGR and, even for b = 1, the most signifi- 
cant deviation is only 2.5-sigma and must certainly not be consid- 
ered as a robust detection of non-Gaussianity. Thus, we find that, in 
the presence of superposed Gaussian signals, the Minkowski func- 
tionals are far less sensitive than the wavelet technique to the pres- 
ence of the non-Gaussian signal due to the cosmic strings map. 

Finally, we perform the Minkowski fun ction als analysis on 
the enhanced noise map discussed in Section 4.4 The results are 
shown in Fig.[l2|. In contrast to the wavelets technique, we find that 
the Minkowski functionals are more sensitive to non-Gaussianity 
in the strings map when the main contaminant is pixel noise as op- 
posed to an EGR. We see that the M2 and M3 functionals provide 
~8-sigma detections of non-Gaussianity at lower limit of the pixel 
temperature distribution. 



6 DISCUSSION AND CONCLUSIONS 

We have investigated the use of wavelet analysis in detecting non- 
Gaussianity in the CMB. Our test images consist of non-Gaussian 
CMB fluctuations from the Kaiser-Stebbins effect due to cosmic 
strings plus some multiple of a Gaussian map with a power spec- 
trum identical to the cosmic strings map. Before being analysed, 
the test images are convolved with a 5-arcmin Gaussian beam and 
Gaussian pixel noise is added with an rms equal to one-tenth that 
of the convolved test map. Using the wavelet technique, we find 
that the non-Gaussian signal can still be detected even when the 
superposed Gaussian map has an rms equal to 5 times that of 
the underlying cosmic strings map. However, statistics based di- 
rectly on the pixel temperature distribution were unable to detect 
the non-Gaussian signal even in the absence of a superposed Gaus- 
sian map. The wavelet technique also produced clear detections of 
non-Gaussianity in the case where the cosmic strings map was con- 
taminated by Gaussian pixel noise with the same rms value; in this 
regime standard methods were once again unsuccessful. 

We also find that the wavelet technique outperforms meth- 
ods based on the calculation of Minkowski functionals. Although 
the Minkowski functional approach yields a significant detection 
of non-Gaussianity for the case in which no Gaussian map is su- 
perposed, we find that, once a Gaussian signal with a equal rms 
value is added, no significant detection is obtained. However, the 
Minkowski functional technique did yield reasonable detections of 
non-Gaussianity in the high pixel noise regime. One possible way 
to improve the results from the Minkowski functional analysis may 
be to introduce into the method some notion of identifying localised 
structure, which is inherent in the wavelet technique. By looking 
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at the distribution of Minkowski functionals as a function of po- 
sition across the map it may be possible to make more significant 
detections of non-Gaussianity. Such a procedure has been investi- 
gated by Novikov, Feldman and Shandarin (1998), who use 'par- 
tial' Minkowski functionals to analyse the COBE 4-year data with 
promising results. 

For both the wavelet and Minkowski functional analyses, we 
have quoted the significance of any detection of non-Gaussianity 
simply as the largest single deviation observed. This is, of course, 
the most conservative assumption that we can make, and corre- 
sponds to the situation in which a high degree of correlation ex- 
ists bewteen the points in the spectrum or between the values 
of the Minkowski functionals at different threshold temperatures. 
Thus, the significances quoted in this paper should properly be con- 
sidered as lower limits. Clearly, if the points were independent, then 
we may simply multiply together the individual significances to ob- 
tain a greatly enhanced detection of non-Gaussianity. In order to 
obtain a value for the true level of significance for the wavelet or 
Minkowski functional approach, it is necessary to derive expres- 
sions for the correlations that exist between the individual points 
in each case. In particular, one would like to define some form of 
likelihood function that includes these correlations. We have not 
carried out such a procedure here, but will investigate these matters 
more fully in a forthcoming paper. 

Finally, we mention how the current wavelet non-Gaussianity 
test might be further improved. The wavelet transform applied here 
decomposes the two-dimensional CMB temperature fluctuations 
into a two-dimensional wavelet basis. The bases used in this pa- 
per were constructed from the direct (tensor) product of existing 
commonly-used one-dimensional wavelet bases (see Section 2.2). 
The main disadvantage of this approach is the mixing of structure 
with different scales j\ and 72 m the horizontal and vertical di- 
rections. As pointed out by Mallat (1989), however, it is in fact 
possible to use one-dimensional wavelet bases to construct two- 
dimensional bases that do not mix scales and can be described in 
terms of a single scale index /. At each scale level j, these bases 
are given by 



The /, wavelet is simply an averaging function at the jlh level, 
while the other three wavelets correspond to structure at the j scale 
level in the horizontal, vertical and diagonal directions in the im- 
age. Since such a basis does not mix structure on different scales, 
we might expect that it would provide wavelet decomposition of 
the CMB map that is more sensitive to the presence of any non- 
Gaussian signal. The application of this form of two-dimensional 
wavelet basis to the detection of non-Gaussianity in the CMB will 
be discussed in a forthcoming paper. 

Even with the above modification, however, the wavelet non- 
Gaussianity test presented here is still restricted to flat two- 
dimensional images. Thus, in the context of the CMB, it can only 
be applied individually to small patches of sky for which the cur- 
vature of the sky is not significant. In order to apply the method to 
all-sky CMB maps, such as the COBE 4-year map or forthcoming 
maps from the MAP and Planck Surveyor missions, it is neces- 
sary to define two-dimensional wavelets on a sphere. A treatment 
of the continuous spherical wavelet transform and its discretiza- 



tion is given by Freeden & Winheuser (1996), and an algorithm 
for calculating bi-orthogonal discrete wavelet bases on a sphere is 
discussed by Schroder & Sweldens (1998). In addition to their pos- 
sible use in detecting non-Gaussianity in all-sky CMB maps, the 
application of such wavelets to the general analysis of CMB data 
may also provide significant advantages in the compression of large 
data sets from the forthcoming MAP and Planck Surveyor satellite 
missions. These topics will also be investigated in a future paper. 
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APPENDIX A: MOMENTS, CUMULANTS AND 
fc-STATISTICS 

For a detailed discussion of moments, cumulants and fc-statistics 
the reader is referred to Struart & Ord (1994) or Kenney & Keeping 
(1954). Here, we outline only the basic definitions. 

Let us consider the probability distribution p(x) of some ran- 
dom variable x. The most straightforward description of p(x) is in 
terms of its moments 

v, = (./), 

where the angle brackets denote an ensemble average. Alterna- 
tively, one could describe the distribution using the central mo- 
ments 

^((x-Vi)'). 

Clearly, if the mean Vi of p(x) is zero, then p,- = v r . The moments 
of a distribution are most conveniently described the its moment 
generating function (MGF), which is given by 



M(h) = (exp(fcc)) 



h 

l+V,/7 + V 2 — + ■ 



•+V r 



h r 



(Al) 



For example, the MGF of a Gaussian distribution with mean Vj and 



variance a is 



M(h) =exp(vih+jG 2 h 2 "< 



(A2) 



Although the moments (or central moments) provide an in- 
tuitively straightforward description of the probability distribution 
p(x), an alternative description in terms of its cumulants has sev- 
eral advantages. The cumulants K r of p(x) are defined through the 
cumulant generating function (CGF) 



K (h) = lnM(/i) 



h z h r 

^l/j + K 2 — + --- + K r — 
2! r\ 



(A3) 



By comparing ( |A1[ ) & (A3), the cumulants of p(x) are given in 
terms of its moments by (Stuart & Ord 1994) 



:(-!)' 



V] 


v 2 


1 


(i) 





1 















Vl 



v 3 




,(A4) 



Evaluating this determinant for the first four cumulants, we find 



K 2 

K 4 



= V! 



v 2 -vf 
v 3 -3v 2 Vi- 
V4 -4V3V1 - 



2V 



- 6v 2 Vi - 3v{ = m — 3po 



There are clear advantages to using cumulants rather than mo- 
ments to describe a probability distribution. Firstly, for two inde- 
pendent random variables x and y, the cumulants are additive, i.e. 



Kr[x + y] ■ 



Secondly, from (A2), we see that the CGF of a Gaussian distribu- 
tion of mean Vi and variance a 2 is simply 

K(h) =v x h+\o 2 h 2 

Thus for a Gaussian distribution K,- = for r > 2. 

From the cumulants and central moments, we can also define 
the dimensionless quantities 



% r = K r /p. r 2 /2 



of which %3 and %4 are better known as the skewness and kurtosis 
of the distribution. For a Gaussian distribution, %3 = and %4 = 3, 
and so it is often useful to define the excess kurtosis X = X4 — 3, 
which equals zero for a Gaussian distribution. 

Let us now consider a sample X\,X2,.-.,Xf/ of size N drawn 
from the distribution p(x). Clearly, these sample values may be 
used to obtain estimates of the cumulants of the parent distribution. 
The most straightforward method is first to estimate the moments 
or central moments of the distribution using 



v r = 



-Yx': 



(=i 

N 



and then use ( |A4| ) to obtain estimates K,- of the cumulants. It is 
straightforward to show, however, that these estimates are biassed, 
so that (K r ) 7^ K r , and this bias is quite pronounced when the sample 
size N is small. It is therefore better to estimate the cumulants of 
a distribution using k-statistics. The first four fc-statistics are given 
by 



K 2 

K4 



Vl 

Nfo/{N-1) 
N 2 fa/(N -1)(N - 



2) 



N 2 [(N + 1)#4 - 3(JV - m\/(N - l)(N - 2)(N - 3) 



The higher-order ^-statistics are listed by Stuart & Ord (1994) up 
to r = 8 but are extremely lengthy to write out in full. Using exact 
sampling theory, it is simple to show that these new estimators of 
the cumulants are unbiassed. 

Finally, we note that, if the parent probability distribution p(x) 
has zero mean, i.e. V] = (as is the case for temperature fluctua- 
tions in the CMB), then there exists an alternative set of unbiassed 
estimators for the cumulants, which are much more straightforward 
to calculate. In this case, the first four estimators are given by 



K 2 
K3 
K4 



Vl 

V 2 



[(N+2)%-3Ml]/(N-l) 
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